#ifndef EQUATIONSOLVER
#define EQUATIONSOLVER

#include "Function.hpp"
#include <iostream>
#include <cmath>
#include <stdexcept> // For std::runtime_error

class EquationSolver {
protected:
    const Function& F; // Reference to the function to solve
public:
    EquationSolver(const Function& F) : F(F) {}
    virtual double solve() = 0; // Pure virtual function
};

class Bisection_Method : public EquationSolver {
private:
    double a, b; // Interval endpoints
    double eps, delta; // Error tolerance
    int Maxiter; // Maximum iterations
public:
    Bisection_Method(const Function& F, double a, double b, 
                     double eps = 1e-7, double delta = 1e-6, int Maxiter = 50) :
                     EquationSolver(F), a(a), b(b), eps(eps), delta(delta), Maxiter(Maxiter) {}

    virtual double solve() {
        if (F(a) * F(b) >= 0) {
            throw std::runtime_error("Function must have different signs at endpoints a and b.");
        }
        double mid;
        for (int iter = 0; iter < Maxiter; ++iter) {
            mid = (a + b) / 2.0; // Compute midpoint
            if (fabs(F(mid)) < eps) { // Check if root is found
                return mid; 
            }
            if (F(a) * F(mid) < 0) { // Root is in [a, mid]
                b = mid; 
            } else { // Root is in [mid, b]
                a = mid; 
            }
            if (fabs(b - a) < delta) { // Check convergence
                if(fabs(F(mid)) > 100 * eps) std::cout << "Warning: May not be the root with function value " << F(mid) << std::endl; //checking if it is a true zero
                return mid; 
            }
        }
        if(fabs(F(mid)) > 100 * eps) std::cout << "Warning: May not be the root with function value " << F(mid) << std::endl;
        return mid; // Return last midpoint if maximum iterations reached
    }
};

class Newton_Method : public EquationSolver {
private:
    double x0; // Initial guess
    double eps; // Error tolerance
    int Maxiter; // Maximum iterations
public:
    Newton_Method(const Function& F, double x0, 
                  double eps = 1e-7, int Maxiter = 8) :
                  EquationSolver(F), x0(x0), Maxiter(Maxiter), eps(eps) {}

    virtual double solve() {
        double x = x0; // Start with initial guess
        for (int iter = 0; iter < Maxiter; ++iter) {
            double fx = F(x);
            if (fabs(fx) < eps) { // Check if root is found
                return x; 
            }
            double dfx = F.derivative(x); // Get derivative value
            x = x - fx / dfx; // Update x using Newton's method
        }
        if(fabs(F(x)) > 100 * eps) std::cout << "Warning: May not be the root with function value" <<F(x) <<std::endl;
        return x; // Return last value if maximum iterations reached
    }
};

class Secant_Method : public EquationSolver {
private:
    double x0, x1; // Two initial guesses
    double eps; // Error tolerance
    int Maxiter; // Maximum iterations
public:
    Secant_Method(const Function &F, double x0, double x1, 
                  double eps = 1e-7, int Maxiter = 50) :
                  EquationSolver(F), x0(x0), x1(x1), eps(eps), Maxiter(Maxiter) {}

    virtual double solve() {
        double x2; // Next approximation
        for (int iter = 0; iter < Maxiter; ++iter) {
            double fx0 = F(x0);
            double fx1 = F(x1);
            x2 = x1 - fx1 * (x1 - x0) / (fx1 - fx0); // Update x2 using the secant method formula
            if (fabs(F(x2)) < eps) { // Check if root is found
                return x2;
            }
            x0 = x1; // Update previous guesses
            x1 = x2;
        }
        if(fabs(F(x2)) > 100 * eps) std::cout << "Warning: May not be the root with function value" << F(x2) << std::endl;
        return x2; // Return last approximation if maximum iterations reached
    }
};

#endif
